************
**Multilevel Models
*************
use "${clean_data}\gfs_cleaned_merged_database.dta", clear
cd "${output}\"

global country_level_controls "v2x_polyarchy zwvs_traditional zwps_index lngdp"
global country_valid "zmort_rate"
global demo "male other_gender secondary tertiary hhinc  age_group2 age_group3 age_group4 age_group5 age_group6 age_group7 foreign_born  employed never_married married live_alone no_children one_child two_children  rural village suburb"
global demo_no_age "male other_gender secondary tertiary hhinc  age foreign_born  employed never_married married live_alone no_children one_child two_children  rural village suburb"
global parent_demo "parents_married living_structure2 living_structure3 living_structure4  child_poverty was_abused"
global extended_controls "parent_relig_index self_relig_index REincome_feelings"

foreach x in $country_level_controls {
gen X_`x'=`x'*parent_relig_index
gen REL_`x'=`x'*PCRQ
}

global interaction_religion "X_v2x_polyarchy X_zwvs_traditional X_zwps_index X_lngdp"
global interaction_pcrq "REL_v2x_polyarchy REL_zwvs_traditional REL_zwps_index REL_lngdp"
gen Xparent_relig_index_PCRQ=parent_relig_index*PCRQ

label var Xparent_relig_index_PCRQ "PCRQ X parental religiosity index"

label var X_v2x_polyarchy "Interaction-parental religiosity X electoral democracy"
label var  X_zwvs_traditional "Interaction-parental religiosity X secular values"
label var  X_zwps_index  "Interaction-parental religiosity X Women's rights"
label var  X_lngdp "Interaction-parental religiosity X GDP per capita"

label var REL_v2x_polyarchy "PCRQ X electoral democracy"
label var  REL_zwvs_traditional "PCRQ X secular values"
label var  REL_zwps_index  "PCRQ X Women's rights"
label var  REL_lngdp "Interact-PCRQ X GDP per capita"

***********
*Finding 1: Predict flourishing--GDP and tradition only country level controls
*by age group
**********

egen age3_cat=cut(age), at(18,30,55,150)
tab age3_cat, gen(age3_group)
label var age3_group1 "age 18-29"
label var age3_group2 "age 29-54"
label var age3_group3 "age 55 and older"

eststo clear

foreach Y in HOPE HAPPY  {
eststo: quietly mixed `Y'  PCRQ $demo_no_age if age3_group1==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo_no_age if age3_group2==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo_no_age if age3_group3==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo_no_age $parent_demo lngdp $extended_controls if age3_group1==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo_no_age $parent_demo lngdp $extended_controls if age3_group2==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  PCRQ $demo_no_age $parent_demo lngdp $extended_controls if age3_group3==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace
}

esttab using "ST1_predict_wellbeing_multi_level_model_by_age.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01 *** 0.001)   nogaps replace

esttab using "ST1_predict_wellbeing_multi_level_model_tstat_by_age.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
b(%15.3f) t(%15.3f) label  star(* 0.05 ** 0.01 *** 0.001)   nogaps replace

***********
*Finding 2: Predict PCRQ--GDP and tradition only country level controls
*by age group
**********

eststo clear

foreach Y in PCRQ  {
eststo: quietly mixed `Y'  $demo_no_age $parent_demo if age3_group1==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  $demo_no_age $parent_demo if age3_group2==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  $demo_no_age $parent_demo if age3_group3==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  $demo_no_age $parent_demo lngdp $extended_controls if age3_group1==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  $demo_no_age $parent_demo lngdp $extended_controls if age3_group2==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

eststo: quietly mixed `Y'  $demo_no_age $parent_demo lngdp $extended_controls if age3_group3==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
estadd local log_likelihood `e(ll)', replace

}

esttab using "ST2_predict_relationship_multi_level_model_by_age.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
b(%15.3f) se(%15.3f) label  star(* 0.05 ** 0.01 *** 0.001)   nogaps replace

esttab using "ST2_predict_relationship_multi_level_model_tstat_by_age.csv", ///
stats(log_likelihood N k, label("Log likelihood" "Sample Size" "k number of parameters")) ///
b(%15.3f) t(%15.3f) label  star(* 0.05 ** 0.01 *** 0.001)   nogaps replace


*******
**Main effect by more detailed age-groups**
********

forval i=1/2 {
gen B`i'=.
gen SE`i'=.
gen UP`i'=.
gen LO`i'=.
gen N`i'=.
}

forval i=1/7 {
mixed HOPE PCRQ  $demo_no_age $parent_demo  if age_group`i'==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
replace N1=e(N) if age_group`i'==1
replace B1=_b[PCRQ] if age_group`i'==1
replace SE1=_se[PCRQ] if age_group`i'==1
replace UP1=B1+(1.96*SE1)
replace LO1=B1-(1.96*SE1)
}

forval i=1/7 {
mixed HOPE PCRQ  $demo_no_age $parent_demo  lngdp $extended_controls if age_group`i'==1 [pw=weight] || iso:, mle nolog vce(cluster iso)  
replace N2=e(N) if age_group`i'==1
replace B2=_b[PCRQ] if age_group`i'==1
replace SE2=_se[PCRQ] if age_group`i'==1
replace UP2=B2+(1.96*SE2)
replace LO2=B2-(1.96*SE2)
}

collapse (mean) N1 B1 SE1 UP1 LO1  N2 B2 SE2 UP2 LO2, by(age_cat)

gen T1=B1/SE1
gen T2=B2/SE2

gen pr1=T1 / (T1^2+(N1-23-1) )^.5
gen pr2=T2 / (T2^2+(N1-28-1) )^.5

drop if age_cat==.

export delimited using "Pooled_PCRQ_Effects_on_Flourishing_by_Age_for_Main_Models_1_3.csv", replace

egen rage=rank(age_cat), track
reg B1 rage
reg B2 rage
reg T1 rage
reg T2 rage
